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We present a method for accurate evaluation of the Green function G(uj,r\, . . . ,rd) at any real 
frequency U) and any lattice vector (n, . . . , rj) for a d-dimensional hypercubic lattice that may have 
anisotropic couplings (Qi, . . . , fid). In this method, we start with an integral representation of G, 
split the oscillatory integrand into combinations of Hankel functions, and deform the integration 
paths into the complex plane to obtain rapidly convergent integrals. We also discuss an alternative 
approach using the Levin collocation method. We report values of the Green function at selected 
frequencies on the branch cut and selected lattice vectors. 



Lattice Green functions arise in numerous problems in 
combinatorics, statistical mechanics, and condensed mat- 
ter physics. In the language of condensed matter physics, 
the Green function G(r, oj) is the propagator for a parti- 
cle of energy (or frequency) ui to move a distance r in a 
tight-binding model on a lattice with a nearest-neighbor 
hopping. Accurate values for lattice Green functions are 
useful as input for further calculations in many-body the- 
ory, such as simulations of Hubbard models^ and non- 
perturbative renormalization group studies P 

Closed-form expressions exist for Green functions on 
many lattices in one, two, and three dimensions, for 
Bethe lattices, and for certain fractal lattices^} see 
Ref. [5]and references therein. In particular, Joyce found 
closed forms for the Green function of the cubic lattice 
at the originp^l a t a general lattice point P and of an 
anisotropic cubic latticeP However, closed forms for hy- 
percubic lattice Green functions in d > 4 dimensi ons a re 
not available in terms of named special functions! 10 * 11 ! 

We recently developed a method^ for obtaining ac- 
curate numerical values for the on-site Green function 
Gq(u>) on a d-dimensional hypercubic lattice. In this 
work we generalize this method to off-diagonal Green 
functions G r (ui) on anisotropic lattices; the generaliza- 
tion is not as trivial as one might expect. We also dis- 
cuss an alternative approach using the Levin collocation 
method to integrate oscillatory functions, and we make 
an interesting observation about the computational com- 
plexity of both methods. Finally, we report values of the 
Green function as test cases against which to benchmark 
implementations or other methods. 

Consider a d-dimensional hypercubic lattice with 
anisotropic nearest-neighbour hopping amplitudes \flk 
(k = 1,2, ...,d) in each direction. The Green func- 
tion G has closed forms in various representations in- 
volving time t, frequency u>, position r, and wavevector 
q = (<?i, . . .,q d ): 



' u> — Qi cos qi — ... — cos qd 
G(q, t) = —iQ(t) exp it(fli cos qi — ... — fid cos q d ) (2) 

G r (t) = i r i+- +r *- 1 e(t)j ri (fi 1 t)...j rd (fi d t). (3) 

The position-frequency Green function G r {u>) can be 



written as a d-dimensional integral over the Brillouin 
zone, 

G ^ = wrL dqi -l dqd 

expifari + . . . + q d r d ) 

ui — fli cos qi — ... — fid cos qd + «0+ 

where we have included an infinitesimal shift in the de- 
nominator as is conventional in condensed matter theory. 
This form makes it clear that G T (ui) has a branch cut 
along the real axis representing a continuous spectrum 
(band) of particle excitations. However, for numerical 
evaluation of G T (u>) it is better to start from the one- 
dimensional integral representatiorP^HIU 

/>oo 

G r (u) =i a dt e lut J ri {flxt) . . . J rd {fl d t) (5) 
Jo 

where a = r± + . . . + r d — 1 and {rk} are integers. 

I. COMPLEX-PLANE METHOD 

The integrand of Eq. ^ is analytic at t = 0, but as 
t — > oo it exhibits irregular oscillations with a slowly 
decaying envelope (t~ d / 2 ), which causes problems with 
traditional quadrature schemes. For example, in d = 4 
the integrand decays as t~ 2 , so it becomes negligible (i.e., 
below machine precision) only beyond t ~ 10 8 . There- 
fore it is useful to adopt the complex-plane technique of 
Ref.lHI 

Split the integral at a point t = T on the real axis by 
writing 

G r (u) = gi +g 2 , (6) 

fji =i a I dt e iut J ri {fl x t) . . . J rd {n d t), (7) 
Jo 

POO 

.92 = i a dte iut J ri {fl 1 t)...J u {fl d t). (8) 

JT 

In Eq. ([8]), apply the identity 

J r {u) = \[H+{u)+H-{u)] (9) 
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FIG. 1: (Color online) Upper panels show integrands for the evaluation of G r (ijj) using the complex-plane integration technique. 
Parameters were d = 4, fife = {1, 1, 1, 1}, cj = 1, ru = {1, 2, 2, 3}. The split point was T = 3. Hue indicates phase and saturation 
indicates modulus. Curves are equal-modulus contours |/| = 1CP 4 , 10~ 3 , . . . , 10 2 . Black circles indicate quadrature points along 
integration path. Lower panels show behavior of each integrand along its integration path; symbols indicate evaluation points. 
It took 165 evaluations of fi, 253 evaluations of /2, and 253 evaluations of f$ to reach 12-digit accuracy using an adaptive 
Gauss-Kronrod quadrature rule. For comparison, the original integral formula Eq. |5| was unable to achieve 6-digit accuracy 
even after 50000 function evaluations. 



where H+(u) = hP(u) and H~(u) 

* are Hankcl functions of the first and second 
kind respectively. Performing a binomial expansion gives 
a sum of 2 d combinations of Hankel functions, so that 



= form 



92 



2 d 



p oo 

J t dte iut l[H%(n k t), 



(10) 



where the sum over {a} runs over all 2 d configurations 
of the Ising variables a k = ±1 and the sum over k runs 
from 1 to d. For \t\ — > oo the Hankel functions behave 
as H°(u) ~ e lau . Thus, for a given {a}, the integrand 
goes as e lwt e *°'i f! i* . . , e »o- d si d t^ anc j j orc | an ' s lemma allows 

us to rotate the integration path into the upper or lower 
half-plane (UHP/LHP) depending on the sign of A[ct] = 
to + oiQ,i + . . . + <Jd^d- Thus 



.92 



2'' 



E 



/ dt e iut Y[H^(n k t) A>0 

•* T k 

dt e iut Y[H%(Sl k t) A<0. 



T-ioo 



(11) 



Together, Eqs. Q, and (11) give a prescription for 
computing G t (oj). 

To appreciate why this procedure is necessary and ef- 
fective, it is instructive to write the Green function in the 



G r (cj) = I dt A(t) +/ dtf 2 (t)+ dtf 3 (t) 

Jt Jt 

(12) 



where 



f 1 (t)=i a e i " t J ri (il 1 t)...J rd (n d t), 



M*) 



2'' 



e n^( n **)' 



2 d 



e n^w, 



(13) 
(14) 

(15) 



and to visualize these three integrations as in Fig. [T] The 
original integrand fi(t) oscillates irregularly, whereas 
/2(t) and /a(r) decay quickly with only a couple of os- 
cillations. 

In order to use standard quadrature routines, it is 
convenient to parametrize the path of integration by a 
real variable r. Returning to Eq. (11) and substituting 
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t = T ±ir leads to 



92 



fdr) 



dr / 4 (t), 



2 d 




+wT II(^r 



where the 2d numbers {%} are 



H±=ff, r ± fc (fi fc (T + *T)), 
= fl£(fi fc (T-tT)). 



A > 
A < 
(16) 

(17) 
(18) 



A naive implementation of Eq. ( 16 1 requires evaluating a 
sum of 2 d products each containing d Hankel functions. 
For an efficient implementation, one should precompute 
the two real exponentials and the 2d complex Hankel 
functions, as described in the pseudocode below: 

Tabulate X± = e ±ur 

Tabulate 7if = H± (ft fe (T + it)) for k = 1, . . . , d 
Set s = 

For each of the 2 d Ising configurations a 

if w + crifti + . . . + <r d n d > o 

Set s = S + X^UkK k 

EXss 

Set s = s-X + Yl k (U^*)* 
End If 
End For 

In order to develop a robust implementation, one 
should consider the following four situations: 

a. Frequency outside band ( Green function not on 
branch cut): If uj > fii + . . . + Cld, then the integrand of 
Eq. ([5| is dominated by the e lu)t term, which decays into 
the UHP. In this case, one can simply rotate the path of 
integration to lie along the positive imaginary axis and 
substitute t = it to obtain 

p oo 

G r H = (-1)''!+-+^ / rf r e - WT / Pl (nir).../ ri (n <l r). 
Jo 

(19) 

This method was put forth in Ref. [13] see also Ref. [9l 

b. Generic frequency within band: In general, the in- 
tegrand fi(r) decays exponentially as t — > oo. This case 
is well handled by standard numerical integration rou- 
tines. 

c. Frequency at van Hove singularity: The Green 
function G r (u>) has singularities at frequencies ±fii ± 
f?2 ± • • • i firf. For a generic anisotropic lattice there 
are 2 d singularities, whereas for an isotropic lattice with 
fife = 1 there are only (d + 1) distinct singularities at 
{— d, —d + 2, . . . ,d}. Whenever w lies at one of these 
van Hove singularities (u) c ), $a{t) decays only as a power 
law, T - d l 2 . This poses no fundamental problem; how- 
ever, canned integration routines that do automatic sin- 
gularity handling may attempt to evaluate the Hankel 



functions at very large u, causing overflow or underflow. 
Therefore, it is advantageous to "fold" the integration do- 
main manually as follows. Split the integral in Eq. ( 16 ) at 
a point r = r\ and apply the transformation r = u /[ ~~ d > 
to the second term to obtain 

92 = f dr U(t) + ^ f ' du u^-^Uin 2 '^). 
Jo Jo 

(20) 

The Jacobian cancels the asymptotic T~ d / 2 behavior, so 
the transformed integrand tends to a constant as u — > 0, 
and can be integrated easily using conventional quadra- 
ture routines. 

d. Frequency near van Hove singularity: If lj f» U) c , 
the integrand fi[r) crosses over from power law decay 
to exponential decay at a large value of the argument, 
t ~ 1/ \cj — cj c \. Any adaptive quadrature routine will be 
forced to sample / 4 (t) in the vicinity of this crossover, 
causing overflow and underflow. In this situation it may 
be necessary to rewrite the integrand in terms of expo- 
nentially scaled Hankel functions e Tlz H^(z), which may 
be computed directly from appropriate recursions or from 
asymptotic expansions such as 



(21) 



There is some flexibility in choosing the split point T. 
With T — oo, one recovers the original formula Eq. d5j). 
Taking T = corresponds to integrating fa{f) from 
to +ioo and fe(t) from to — ioo respectively. This is 
an attractive choice because Hankel functions of imagi- 
nary argument can be computed cheaply by decomposing 
them into Bessel K and Bessel / functions.^ However, al- 
though the original integrand f\{t) was analytic at t = 0, 
the split integrands /2(f) and fz{t) are singular at t = 0. 
In the case of the on-site Green function (r = 0), /2(f) 
and have logarithmic branch points at t = 0, which 
are weak, integrable singularities, so it is permissible to 
rotate the integration paths as was done in Ref. 12, How- 
ever, for r 0, f2(t) and fs(t) both have a pole at t = 
(which may be a high-order pole), so that it is illegal to 
rotate the path about t = 0. Then it is necessary to use 
a split point T > sufficiently far away from the pole at 
t = 0. 

The optimal value of T depends on the implemen- 
tation of Bessel and Hankel functions and the quadra- 
ture routines available in one's software package. In 
the current version of Mathematica at the time of writ- 
ing (8.0.1.0), the special function routines are slow 
(HankelHl [0 , 2 . +3 . 1] takes almost 0.7 milliseconds) and 
inaccurate (the error in BesselJ [0 , 11 . ] is 1000 times 
larger than machine precision), so we have not attempted 
detailed fine-tuning. 

There is also some flexibility in the choice of integration 
paths. Instead of integrating fc(t) along the straight line 
T — > T + ioo, it may be desirable to use the path T — > 
iT — > ioo, because the integrand is cheaper to evaluate 
when t is pure imaginary. 
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For general r, the binomial expansion leads to 2 d terms. 
In certain cases these are not all distinct. For example, 
for r = on an isotropic lattice, the binomial expansion 
of Jo(t) d only gives rise to (d + 1) distinct terms^. 



II. LEVIN COLLOCATION METHOD 

As remarked earlier, the integrand of Eq. ^ has irreg- 
ular oscillations that make it difficult to integrate using 
traditional quadrature schemes. However, integrands of 
this form are amenable to relatively modern integ ration 
techniques such as the Levin collocation method^^. 
Mathematica 8.0.1.0 is able to perform an automatic 
symbolic analysis of the integrand to determine a suitable 
set of basis functions and carry out the Levin method. 
We summarize the method below, assuming the summa- 
tion convention for brevity. 

The Levin method deals with an integral of the form 
/ = J dt fiWi (summed over i = 1,2,3, ...,n) where 
fi(t) are slowly varying envelope functions and Wi(t) 
are oscillatory basis functions. The n basis functions 
satisfy a matrix ordinary differential equation (ODE), 
w[ = AijWj, where Aij(t) vary slowly. The fundamen- 
tal theorem of calculus gives the integral in terms of the 
antiderivative of the integrand: 



I = Fi(b)wi(b) - Fi(a) Wi (a) 



(22) 



where (i^iOj)' = fiW t . Simple manipulations show that 
the antiderivative envelope functions Fi (t) satisfy the ma- 
trix ODE 



F[ + A 3l F 3 



; 'kernel" A^t) 



= fi, (23) 
and "forcing" fi{t) vary 



in which the 
slowly. 

The general solution for Fi(t) contains an admixture 
of up to n oscillatory "modes" . However, Levin ob- 
served that there is one solution that is slowly varying, 
and that this "particular solution" can be well approx- 
imated by solving the ODE using a collocation method 
as follows. Choose a set of collocation points {ti} for 
I = 0, 1, 2, . . . , to — 1. It is convenient to use the ex- 
trema of the Chebyshev polynomial of degree to (suitably 
scaled): U = |[6+ a — (b — a) cos ^zrf]- Choose an ansatz 
for Fi(t). It is convenient to use a linear combination of 
Chebyshev polynomials of degree less than to, weighted 
by coefficients Ci k (k = 0, 1, 2, . . . , to — 1), 

(24) 



Fi(t) ~ c ik u k (t) where u k (t) = T k ( 



2t-b- 
b—a 



The derivatives are series in Chebyshev polynomials of 
the second kind, U k : 



F[{t) a c lk u' k (t) where u' k (t) = ^C7 fc _ x (2^). (25) 
Let us demand that the ansatz satisfy the ODE exactly 



at the collocation points. Putting Eqs. (24 1 and (25 1 into 



Eq. (|23j) leads to 



Aji(ti)u k (ti)] c jk 



Ml). 



(26) 



This is a system of mn linear equations, which can be 
solved to give the mn coefficients Cj k in 0(m 3 n 3 ) time. 
From the numerical values of the Cj k one can then use 



Eq. ( 24 ) in Eq. ( 22 ) to estimate the integral /. 

Let us now examine the application of the Levin 
method to the computation of the Green function 
G rir2 „. rd (ui) via Eq. (|5|. In general, this requires 2 d ba- 
sis functions, Wi(t), which are all possible products of 
Bessels and Bessel derivatives, 



J ri {Slxt)J r2 {Q, 2 t) ■ ■ 
J' ri {nxt)J r2 {9. 2 t).. 
J ri (fiit)J' (fi 2 t).. 



J rd (Q,dt), 
J rd (Qdt), 
J rd (p,dt), 



4 i (n 1 t)j^(n 2 t)...4 d (n d t). 



(27) 



For clarity of exposition we illustrate the procedure for 
the on-site Green function on an isotropic lattice in d = 4, 
Gq(uj) — —i J e lult Jo(t) 4 , which requires only n = d + 
1 basis functions. The basis, forcings, and differential 
matrix are 



b»i} = 



V e-*Jo(t) 4 ) 



{fi} = 







W 





-4 








o\ 


i 


iw-f 


-3 











2 


IU1 — \ 


-2 











3 




-1 


V 








4 


ioj J 



(28) 



Using Eqs. (|24|, ([25} , (|26j), and ([22]) with to = 11 colloca- 
tion points gives the integral from t — 10 to t = 100 with 
absolute error about 10 -6 (relative error about 10 -3 ). 
The intermediate quantities are illustrated in Fig. [2j 

It is interesting to note that in the Levin method, the 
properties of the integrand are encapsulated in the A$j (t) , 
so that the integrand is never evaluated at any t (unlike 
in traditional quadrature methods). The special func- 
tions e lwt , Jo(i), and J' {t) are only evaluated at t = 10 
and at t — 100. Thus, the computational time is dic- 
tated not by the number of integrand evaluations, but 
by the solution of the mn simultaneous equations for the 
Chebyshev coefficients. 

The above-described Levin method has some draw- 
backs in the present context: 

• Integrating from t = to t — oo requires some 
extensions due to the singularity in A^jit —¥ 0) and 
the infinite upper limit. 

• Solution of the collocation equations may be less 
controlled (with respect to roundoff accumulation) 
than performing quadrature on the smoothly de- 
caying integrand in Eq. (16). 
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FIG. 2: (Color online) Graphical examination of the Levin method as applied to the integral I = dt e 1,5li Jo(i) 4 . The 
oscillatory integrand is written as a combination of n — 5 basis functions Wi(t) modulated by integrand envelope functions 
fi(t) (both real and imaginary parts are plotted). The antiderivative (indefinite integral) is represented in the same basis by 
antiderivative envelope functions Fi(t). The Fi(i) are approximated by Chebyshev polynomials and the Chebyshev coefficients 
Cjk are determined by demanding that Fi(t) satisfies a matrix ODE at m — 11 collocation points t;. The values of fi(ti) and 
Fi(ti) are illustrated as circles (although the latter are not needed for the calculation of 7). The change in the antiderivative, 
-Fi(100)u>i(I00) — Fi(10)wi(10), gives a good approximation to I. (We are using the summation convention.) 



• Since the true Fi(t) is smooth, the simplest way 
to improve collocation accuracy is to increase the 
number m of Chebyshev polynomials, but this is 
restricted by the 0(m 3 n 3 ) time complexity for 
solving the collocation equations. In contrast, 
Clenshaw-Curtis quadrature for the smoothly de- 
caying integrand in Eq. (16), which essentially in- 



volves integrating the Chebyshev interpolant, can 
be taken to arbitrarily high Chebyshev degree 
m, since the complexity scales only as 0(m) (or 
O(mlnm) if the quadrature weights are computed 
on the fly). 

• When lo lies at a van Hove singularity, the Levin 
method experiences problems due to resonance be- 
tween the e Iwt and Jo(t) factors: the integrand ac- 
quires a non-oscillatory component that needs spe- 
cial treatment. 

We observe an interesting parallel. The complex-plane 
method for Go(w) involves a sum of d + 1 terms, and for 
G ri ... rd (uj) the sum involves 2 d terms. The Levin method 
for Go(w) requires n = d + 1 basis functions, and for 
G ri ... rd (uj) it requires n — 2 d basis functions. However, 
the complex method involves quadrature, which takes 
time linear in d+ 1 (or 2 d ), whereas solution of the Levin 
collocation equations takes time cubic in d + 1 (or 2 d ). 
Thus for large d the Levin method is more expensive. 

Taking all the above considerations into account, we 
feel that the complex-plane method is more suitable for 



evaluating hypercubic lattice Green functions of the form 
Eq. ([5]). Nevertheless, one's choice may be influenced by 
the availability of complex Hankel functions and quadra- 
ture routines in one's preferred language or software plat- 
form. We have not attempted an exhaustive study of 
whether other methods for oscillatory integration^!^ 
are able to deal with Eq. ([5| . 



III. SELECTED VALUES OF G 



The following values of the Green function for isotropic 
cubic and hypercubic lattices, i.e., with {f^} = (1,1,1) 
and (1,1,1,1) respectively, were computed using the 
complex-plane method implemented in Mathematica (see 



() 



auxiliary file HLGFPackage . nb) : 



Gooo(3) 
Gioo(3) 
G U0 (3) 
Giu(3) 
G2oo(3) 
Gooo(O) 
Gioo(O) 
Gno(O) 
G m (0) 
G 200 (0) 



0.50546201972 

-0.17212868638 

0.11038286738 

-0.08715670880 

0.08577862908 

-0.89644078878i 

0.33333333333 

0.18578752146; 

-0.27566444771 

0.15329070292^, 



(29) 



As expected, ImG r (ti;) = when \uj\ > d (the 
density of states is zero outside the band). By 
symmetry, ReGn(0) = 0. The Green function 
obeys the discrete Hclmholtz equation, uG r {uj) + 

\ Y?k=i X^=±i G r+aek {io) = S r , where {e k } are unit 
vectors and S T = 1 if r = and zero otherwise. At 
the band bottom (w = — d) the Green function obeys 
the discrete Laplace equation, so ImG r (- d) = and 
ReG r (— d) < Vr. The Green function at the top of the 
band alternates in sign between adjacent sites. 

The accuracy of the above results (11-13 digits) was 
limited partly by the large errors in the Bessel-related 
functions supplied by Mathematica 8.0.1.0. 



G 00 oo(4) = 0.309866780462 

Giooo(4) = -0.05986678046 

Gnoo(4) = 0.02542940754 

Gmo(4) = -0.01546809528 

G mi (4) = 0.01118185767 

G 20 oo(4) = 0.01649101798 

Goooo(0) = -0.90272857832; 

Giooo(O) = 0.25000000000 

Gnoo(O) = 0.15098515279; 

Gmo(O) = -0.10132118364 

G im (0) = -0.20025275758; 

G 2 ooo(0) = -0.00318233840; 

G ooo(l) = 0.3726972107993 - 0.6681496264378; 

G 0000 (2) = 0.5680714850367 - 0.3573566432144; 

G 0000 (3) = 0.4358824699995 - 0.1063899831047;. (30) 



IV. CONCLUSIONS 



We have presented an efficient method 
[Eqs. ||6[),(|7]),(|l6"))], based on complex analysis, for 
calculating Green functions of d-dimensional anisotropic 
hypercubic lattices at arbitrary lattice vectors. We also 
discuss the merits of an alternative approach using the 
Levin collocation method, but we provide reasons to 
prefer the complex-variable method. 
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